Quasi-stationary Stefan problem and computer simulation of 

interface dynamics. 
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Abstract. The computer simulation of quasistationary Stefan problem has 
been realized. Different representations of the Laplacian growth model are 
considered. The main attention has been paid for the interface dynamics 
represented by integro differential equations. Numerical approach has been 
realized by use of interpolating polynomials and exact quadrature formulae. 
As a result system of ordinary differential equations has been simulated. 



There are various physical processes in nature, such as crystallization, combus- 
tion, deposition etc., which make it is possible to single out moving boundary having 
physical properties different on its either side. The location of the given boundary is 
not known beforehand and is determined self-consistently by the distribution of its 
physical parameters in the ambient space. Such problems as usually are singled out 
as a separate section of modern mathematical physics, where the problems with the 
moving boundaries, also known as Stefan problems, are considered. The number of 
publications, devoted to this problem, is rather great. Therefore, we will mention 
only some monographs H, Q including a lot of bibliographies on mathematical 
models of various physical systems containing moving boundaries. 

Such problems, as usual, have no analytical solution due to the essential non- 
linearity stipulated by representation of boundary conditions for selfconsistent pa- 
rameters allocated on the free boundary. Though it is sometimes possible to find an 
analytically self-similar solution to a Stefan problem in one-dimensional or quasi- 
one-dimensional case, concerning the two-dimensional case such problems can be 
solved using computer simulation only. Due to the fact that problems with the mov- 
ing boundaries describe a wide class of important technological processes, such as 
crystallization, melting, etching, oxidation, deposition, diffusion etc., development 
of numerical methods of their solution is especially important [[j], |2| . 

1. Model 

Let us consider crystallization of an undercooled melt located in infinite con- 
tainer, with initial temperature T < T c , where T c is the temperature of crystalliza- 
tion H . It is supposed that growth of the interface of the solid and liquid phases T 
is controlled by the extract of heat generated in the process of crystallization inside 
the melt (Fig. |l|). Then the velocity v n of the motion of the crystallization front 
along the normal vector n to the interface T is given by the expression: 
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Figure 1. A schematic view of the interface during crystallization. 



(1.1) v n = -JI(vT| r n). 

In the melt the distribution of temperature T is determined by the heat equation: 

dT 

(1.2) c pP — = xAT. 

Here c p and p are the specific heat and density of the liquid phase; xr is thermal 
conductivity; L is the latent heat of crystallization. The temperature T s at the 
crystallization front usually satisfies the Gibbs-Thomson condition: 

(1.3) T S =T C (1-Kd ), 

where K is the local curvature of the surface; do is the capillary length proportional 
to the surface energy. 

It is known that in systems described by equations ( [1.1] ) - (1.3) regular dynam- 



ics of the interface is unstable. The basic issues for understanding the processes 
of dissipative structures formation during crystallization are presented in 0, 0, 
clarifying the reason of formation of these structures for instance the morpholog- 
ical instability. In case of the undercooled melt the mechanism of this instability 
is stipulated by magnification of the motion velocity of convex segments of the 
surface penetrating into the undercooled region. The latter stimulates their faster 
propagating into depth of the melt. 

The realization of instability of the homogeneous state of an interface for crys- 
tallization of the undercooled melt, within a wide interval of wave numbers and the 
essential instability of the whole process of shaping, reduces itself to large varieties 
of possible structures and next to the reproduction of the self-similar forms of the 
interface in smaller scales. It also indicates the fractal character of their formation. 
In spite of apparent simplicity of the problem statement, even in two-dimensional 
domains it is often senseless to perform direct numerical calculations for solving the 
problem. 

In same cases the microscopic temporal scales of the parameters specifying the 
problem changes can appear to be considerably smaller than the temporal scales 
of dissipation and distribution of an appropriate unlocal field. As a result the 
distribution of the unlocal field can be considered quasistationary subject to the 
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location of the boundary. In this case one can construct the equations of surface 
motion containing only undefined fields specified on the considered surface. It is 
the basic methodology for this approach. 

The essence of the method devised in |6| implies replacement of a consistent 
system of equations of motion by means of more simple ones, while keeping all 
the properties of the origin problem. As the velocity in the obtained equations is 
defined locally at each point, we obtain a set of equations of local dynamics. Some 
simplifications, however, restricted application of the approach proposed. We can 
see that the equations obtained in Q are satisfied for undercooling, when the 
characteristic scale of an unlocal field (for example temperature) It = x (cpv) 1 
is much smaller than the radius of surface curvature, r i.e.: l t <C r. Although in 
this case one can not obtain solutions without numerical simulation making use of 
appropriate approximations allows to find new solutions, for example those in the 
form similar to snowflake. 

If It r , then in the field of spatial changes of the temperature it is possible 
to allocate two characteristic subregions. The first one, with characteristic scale, 
comparable to l t , corresponds to the change of the temperature, due to the effect of 
thermal conduction, when the curvature of the boundary element can be neglected. 
The second subregion lies in proximity to the boundary element, and its thickness 
is equal to the average curvature radius r. The change of the temperature and its 
gradient in this subregion depends on curvature of the boundary. The modulation 
of the surface temperature and its spatial distribution, as the first approximation 
in the parameter r/l t , corresponds to the Laplace equation. Since harmonic func- 
tions as the solution of Laplace equation describing the distribution of temperature, 
methods of complex analysis appear to be very effective. This book presents meth- 
ods of conformal mappings for analysis of motion of boundary surface at least in 
two-dimensional case. Thus, if It 3> r the problem of crystallization of undercooled 
melt ([O]) - fll-3| ) and the problems considered above are formally equivalent and 
are typical representatives of the quasistationary Stefan problem. Making use of 
methods of complex variables allows us to reduce the problem discussed to a system 
of ordinary differential equations ]8|, |9], |l(| . 



In this case problem (1.1) - ( ]1 ,3| ) becomes one of the simplest among the ones 
considered above and is a kind of aquasistationary Stefan problem. However, even 
in this simple case, determination of the boundary is a rather complicated mathe- 
matical problem. 



Reformulate now the problem (1.1) - (1.3) in terms of Laplacian growth prob- 
lem. We investigate the growth of the two-dimensional region Q and suppose that 
its boundary is T, which represents the physical interface (Fig. |l|). The field <p 
outside the region Q, satisfies the Laplace equation 

(1.4) A<^ = 0. 

The boundary T grows at a rate, that is proportional to the normal gradient of the 
field of the interface. Therefore the evolution of the interface is governed by the 
following equation 

(1.5) v n = -V n </j| r . 

The potential field tp satisfies the following boundary condition on the interface 

(1.6) cp\ r = d k\ r , 
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where k is the curvature of the interface T, and do is the dimensionless surface 
tension parameter. 

For this type of system the utilization of the theory of function of complex 
variable is very fruitful. Theoretical methods of the functions of complex variable 
appear to be rather useful in the construction of the effective surface dynamics of 
the interface at least in the two-dimensional problem. In this case the central issue 
of this problem belongs to the conformal mappings of the physical regions of the 
space onto some special kind of regions. A large number of models for interfacial 
dynamics had been proposed up to now. Reviews of these problems can be found in 
11, |12fl . The most famous equation was introduced by Shraiman and Bensimon 



14 1 . It is widely known as the conformal mapping equation and is still under 
investigation. 

In examples of the quasistationary Stefan problem, the main point of its solution 
lies in construction of an appropriate potential specified by its boundary value. 
Physically it is possible to do only when the boundary of the domain of the potential 
is rather simple. The simplest examples of such region is a disk. If we know a 
conformal mapping of a disk on some region, the solution to the posed boundary 
problem reduces to change of variables in the solution for the disk. Thus, solution 
of the quasistationary Stefan problem for the two-dimensional region is practically 
reduced to determination of an appropriate conformal mapping or its equation 
of motion. It is possible to compare the non-stationary Stefan problem with an 
equivalent quasistationary Stefan problem. It allows us to solve a wider class of 
problems in framework of the equations of motion of the appropriate conformal 
mappings. 

Let us consider the approach developed by Shraiman and Bensimon jj^, 12, 
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This approach is based on the Riemann mapping principle. This 
principle states the existence of a conformal map from the exterior of domain Q 
onto the standard domain, for example, the exterior of unit disk, with the boundary 
r corresponding to unit circle. In the outline under regard one can parametrize 
evolution of the interface in z plane by means of time-dependent conformal map 
F (w, t) which takes the exterior of the unit disk at each moment of time, \w\ > 1, 
onto the exterior of Q. Therefore the evolution of the interface can be presented as 
follows: 

(1.7) r(0,t)= lim F (w, t) , < 6 < 2tt. 



As it is shown by Shraiman and Bensimon [14|, in the case of free surface tension 



(do = 0) the boundary map satisfies the following equation of motion 



dT(9,t) _ . dT(9,t) 

(1 ' 8j ~d^~ l ~d9~ b 



dT(9,t) ~ 2 



09 



where S [■] means the Schwartz operator |13|. For the exterior of a unit circle with 



some boundary condition / (9) the Schwartz operator can be presented as follows: 

(1-9) S [/ (9)} = - — / dQ'-^-f (9') + iC, 

2ir e w — w 



where C is an arbitrary constant. 
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It was shown p. 5| that equation (1.8) is equivalent to equation of motion of 



certain field r\ (0,t) 



ar(8,t) 



ae 



dn{9,t) _ 1 dr)(6,t) 1 



dt R[t) 89 2tt 



d8'v„(6',t)r)(6',t) cot ■ 



R(t) 



v(0,t)v n (e,t)-— / ^e' \t)v n (e' \t)dff' 



+ 



r) 2 (6,t)v n (e,t) d i 

R(t) 89 2ir 
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dO' lnr)(0',t) cot ■ 



(1.10) 



V(0,t) d 1 
R(i) d6 2Ti 



d6'v n {6',t)r](6',t) cot ■ 



considered in [12, 16, 17]. Here the radius R and normal velocity v n satisfy to 
the next equations: 

(1.11) ™ = ± 

V ' dt 2tt 



2tt 



d9'n(6',t)v n (9',t). 



d6' cot x 



x^-v(6',t)do(9',r)) w 



— I d9"v<>\ ■ 

2tt 



o 



In 77(0", t) 



+ 



(1.12) 



+ 



^ s d 1 
U(t) 96» 2vr 



d e' cot —-r ] (e',t)d (e',r ] ) 



In local approximation last relationship has the form |12|, |l6], 17 1: 



(1.13) 



v n (9,rj,t) 



R(t) \ R(t)86 2 



The equations ( 1.10 ), ( 1.11 ), ( 1.13 ) are one-dimensional equations determined 
on the interval [0, 2tt]. 

Now let us consider the linear stage of stability of a cylindrical boundary. 

In this case it is natural to put dp (9, rf) = d a = const, /3 (6,rj) — f3 = const 

. Linearizing equations (1.10), (1.11), (1.13) with respect to perturbations Stj — 
t 

5f] Q exp(J 7 (t') dt' + ikd) at k = ±1, ±2, ... and taking into account the relation 




2 71 



i / dd' cot • 

2tt 



jke _ ■ Ake 



i e sign/c 
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we obtain the following characteristic equation 



(1.14) 



7 




R 



fc_2 + ^fc 2 (l-fc) 
R 



It can be seen from equation (1.14) that the instability is realized for rather 
short-wave fluctuations with k > 2 IE]. At d ^ the fluctuations with wave 
number 




(1.15) 



have the maximum increment of increase. The instability as do ~~ > is observed 
within a wide spectrum of wave vectors, which causes complicated fractal structure 
of dendritic formations. The surface energy damps the instability subject to chosen 
wave vectors and gives rise to the minimum radius of curvature of a dendrite. 

2. Numerical simulation of the interface dynamics 

In the present section numerical simulation of the surface dynamics is given. 
We consider the properties of the interface dynamics from the viewpoint of singular 
integral equations. 



2.1. Properties of the interface dynamics. Equations (1.10) and (1.11) 
are not directly related to the physical mechanism determining motion of the in- 
terface in real space. They determine the variation law governing the mapping of a 
circle into the contour T under the given normal velocity v n . In framework of this 



approach we consider properties of equations (1.10), (1.11). 

We will investigate the periodical structure formation, supposing that the struc- 
ture will have period 27r/fc, where k is an integer. This make it possible to consider 
the solution only on the interval [0, 27r/fc]. In fact, we shall consider some integrable 
periodic function 



(2.1) 



f(9 + n27r/k) = f(9) 



In this case the any singular integrals J in the equation ( 1.10 ) will look like: 

2ir/fc 



2tt 



1 



(2.2)— / dO'f (6') cot — 



(2.3) 
where 



n=0 - 



2ir/k 



d6'J(e')U{9' - Q)dv! 



k-l 



»(u' -«) = £ cot t ± n2 l' k 



So, we can use for computer simulation the interval [0,2-7r/fc]. 
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2.2. Numerical solution of the interface dynamics. The derived sys- 
tems of nonlinear integro-differential equations were investigated numerically. The 
numeric simulation was based on approximation of unknown functions by trigono- 
metrical polynomials. Expansions of this type make it possible to use well-known 
quadrature formulae to evaluate singular integrals. As a result, the system of equa- 
tions (1.10), (1.11) and (1.13) reduces to a system of ordinary non-linear differential 
equations with filled Jacobian. 

During computing experiments the approach based on approximations of searched 
functions by trigonometrical polynomials 



18 



IS 



(2.4) 



fn(u) 



^ 2n-l 

r E f(uj)sm[n( 



2tt 



u — u, cot • 



3=0 



with the given nodes set 



(2.5) 



*3 



j=0,l,...,2n-l. 



was used. The representation of the function / (u) in form (2A) makes it possible 
to apply the obtained analytical expressions practically for all derivatives of the 
function / during calculations. 

Besides, expansion of the function / (it) into a trigonometrical polynomial 
makes it possible to use quadrature formulae | 
integrals with Hilbert kernels, that is 



19 1 for calculation of singular 



2n-l 



(2.6) 

3=0 

being defined on the system of nodes 

2nm + 1 



cot ■ 



(2.7) 



2?r 



m = 0, 1, 2n — 1. 



Taking now in to ac count that all integrand functions on the right hand side of 
equation s (px| ), ( |t| ) are of the same kind it is possible to present them in the 
form of (2J3) making the following assignments 



(2.8) 



m = 0, 1, In — 1. 



As a result from the equations of interface dynamics we obtain the system of ordi- 
nary differential equations: 



(2.9) 



dt 



2n-l 



cot ■ 



where F m is a nonlinear operator being determined by the right hand side of equa- 
tions ( 1.10 ), ( 1.11 ) ( 1.12 ), fj can take values v n (uj)r](uj),\nr](uj),r](uj) and etc. 

Numerical integration of specified systems was carried by means of discrete 
calculation of the Jacobian when realizing implicit methods of Girr and Adams 
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Figure 2. Non-uniform distribution of the r](6,t),ip(9,t) and 
the corresponding interfaces for the following initial perturbation: 
r)(9, to = 0) = 1.0 — 0.01 cos(46>); do = 0.05; interfaces on the figure 
correspond to the time steps (t ste p = 0.5 *10 4 ) 18 times are plotted, 
spaced between < t ^ 0.8 * 10 4 . 



|20j| . Choice of the value of discretization and the local error of integration was 
made so that the necessary accuracy of the searched solution and the stability of 
calculations were ensured. Making use of variable r\ changing in time, subject 
to the transformations 

dy = *@* f 6*(o + m + 
n[6,t) \ 



(2 . 10) * . ^^ +m+ i) 



for the system of equations (1.10), ( ]1 . 1 1 ) and (1.13) it is possible to reproduce 



geometrically defined interfaces in parametric form. 

In this text we evidently do not claim the complete investigation of pattern 
formation in free boundary problems. Our goal is to show that developed approach 
can be applied to investigation of interface evolution. Results of computer simula- 
tion presented below give rise only to some sketch of the structures realizing in the 
systems under consideration. 

During numerical simulation evolution of the interface was investigated. The 
product term j3 affect only on time scale of the process. Therefore, in the analysis 
under regard (3 = 1 is taken. It was found out that the surface energy influences the 
behavior of the field rj(9,t) and interfaces too much. At small values of the surface 
energy rather small-sized structures can be formed. As do increases the structures 
take a more smooth form. 

Figure |^ illustrates the characteristic form of non-uniform distributions of the 
field r\ (plot a) and the corresponding interfaces (plot b). 

It is possible to observe that relatively smooth variations of the interface cor- 
respond to rather great variations in the field rj as well as in the angle ip. Thus, 
relatively small errors in determining -q do not cause great changes in the interface 
geometry. In addition, as is seen from Figures || relatively even long sections of 
the interface shrink when mapped onto the interval [0, 2w], whereas sharp changes 
in the contour T, conversely, are stretched out. In this sense such mappings are 
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Figure 3. Comparison of interface profiles for different values of 
surface tention do : do = 0.15— a); do = 0.12— b); do = 0.09— c); 
do = 0.06-d); d a = 0.04-e); d = 0.03-f). Initial pertur- 
bations for all the figures are the same: r)(6,to = 0) = 1.0 — 
0.01cos(4#); (t step — 0.5 * 10 4 ) 18 times are plotted, spaced be- 
tween < t < 0.8 * 10 4 . 



"adaptive" with respect to the information about the detailed geometric structure 
of the interface. 

The characteristic geometry of the structures obtained as a result of computer 




The smaller 
are different 



simulation of a small perturbation is presented on Figures 
capillary length is the more diverse structure arises. Figures 
from each other due to periodicity of the initial perturbation. 

If the system is influenced by an uncorrelated noise with the normal law of 
distribution unregular structures arise. During simulation it was taken that on 
each step of integration the velocity at point i of the interval [0, 27r] equal to 

1 /2 

v n (0i) + (t)/R) Sv n (0i), within 5v n (9i) being chosen in the correspondence with 
the normal law of distribution in limits (— 10~ 5 , 10~ 5 ) (Fig. ||). 

For simulation of anisotropic growth of crystals the following expression (similar 
as in 



21 1) was used: 



d (0) = (1 - 157 4 cos 47)d = 



(2.11) 



271 



1 — 157 4 cos 4 



1 

27 



d6»'ln?7(6»')cot' 



The interfaces corresponding to the capillary length (2.11), subject to the 
anisotropy of surface energy have the form of truly looking cross (Fig. |6[ and 
0), despite that at the initial moment the boundary surface presents a circle (up to 
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Figure 4. Comparison of interface profiles for different values of 
surface tention do: do = 0.15 — a); do — 0.12-b); do = 0.09— c); 
do = 0.06— d); do = 0.04— e); do = 0.03— f). Initial perturbation for 
all the figures is the same: rj(6,to = 0) = 1.0 — 0.01 cos(6#); (t s tep = 
0.5 * 10 4 ) 18 times are plotted, spaced between < t ^ 0.8 * 10 4 . 



Figure 5. Interface dynamics stimulated by external noise. 



8v n \t=o — 10~ 5 ). According to the value of 7 4 such an anisotropy is exhibited at 
small or large values of R and t. The realized distributions little depend on initial 
perturbations and reduce frequently to rather complicated and exotic surfaces. A 
more developed structure in the case of anisotropy is presented by Fig. §1 An 
example of structures arising from nonhomogeneous initial conditions is illustrated 
by Fig. [| a). In the same time Fig. ^, b) illustrates an example of structures 
developing from a flat boundary. 

It should be noted that at do — the results concerning smooth boundary 
under stochastic excitation of short-wave perturbations, caused for example by 
errors of calculations do not lead the developed structures due to instabilities of 
computation. 
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